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Crossover from Positive to Negative Interlayer Magnetoresistance in Multilayer 
Massless Dirac Fermion System with Non- Vertical Interlayer Tunneling 

Takao Morinari* and Takami Tohyama 
Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan 

We present a theoretical description of the interlayer magnetoresistance in the layered 
Dirac fermion system with the application to the organic conductor a-(BEDT-TTF)2l3 un- 
der pressure. Assuming a non- vertical interlayer tunneling and including higher Landau level 
effects we calculate the interlayer conductivity using the Kubo formula. We propose a physical 
picture of the experimentally observed crossover from the negative interlayer magnetoresis- 
tance, where the Dirac fermion zero-energy Landau level plays a central role, to the positive 
interlayer magnetoresistance that is the consequence of the Landau level mixing effect upon 
non-vertical interlayer hopping. The crossover magnetic field depends on the Landau level 
broadening factor and can be used to determine the Dirac fermion Landau level energy 
spectrum. 
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1. Introduction 

In condensed matter systems, it is possible to realize a massless Dirac fermion like excitation spectrum 
of electrons with the speed of light being replaced by a Fermi velocity. For example the electron energy 
dispersion described by the tight-binding model with nearest neighbor hopping on the honeycomb lattice 
are well approximated by a massless Dirac fermion spectrum at the high-symmetry points at the corners of 
the first Brillouin zone. 1 ) Since the discovery of unconventional integer quantum Hall effect in graphene, 2,3 ) 
which is a two-dimensional carbon material with a honeycomb lattice, Dirac fermions realized in condensed 
matter systems have attracted much attention. 

In conventional two-dimensional electron systems at low temperatures and strong perpendicular mag- 
netic fields, the Hall conductivity a xy exhibits plateaus at ne 2 /h with n the number of filled Landau levels. 4 ) 
By contrast in two-dimensional Dirac fermion systems a xy exhibits plateaus at a xy = ±(\n\ + l/2)e 2 /h 
with n the Landau level index. 2 ' 3 ) (Here we do not include the degeneracy factor due to spin and valley.) 
In the Dirac fermion system the energies of the Landau level are described by 



E n = sgn(n)Cv^p, (1) 

with n = 0, ±1, ±2, ... and B being the magnetic field. The shift of e 2 / (2h) in the Hall conductivity plateaus 
is associated with the existence of the zero energy Landau level, n = 0. 5 ^ 
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Another massless Dirac fermion system is realized in a-(BEDT-TTF)2i3 under pressure. 6-8 ) (Here 
BEDT-TTF denotes bis(ethylenedithio)-tetrathiafulvalence.) This system has a layered structure, and the 
conducting layers consist of BEDT-TTF molecules. Between conducting layers there is an insulating iodine 
layer. The system exhibits small inplane conductivity anisotropy while the conductivity perpendicular to 
the plane is much less than that in the plane. 9 ) This strong anisotropy suggests that the system consists 
of weakly interacting conducting layers. In the ambient pressure, a metal-insulator transition occurs at 
135K 9 ) associated with a charge ordering. 10 ~ 13 ) This metal-insulator transition is suppressed by applying 
pressure. 14,15 -* Under pressure larger than ~ 2GPa, the system is metallic and the resistivity is almost 
temperature independent while the Hall coefficient shows strong temperature dependences. 16 ) Kobayashi et 
al. calculated the energy dispersion 6 ' 7 ) using the tight-binding model. The transfer integrals were estimated 
from X-ray diffraction experiments. 17 ) They found that the electronic band structure near the Fermi level 
is described by a tilted and anisotropic Dirac cone. Such a Dirac cone structure was also obtained by first 
principles calculations. 18 ' 19 ) 

The presence of the Dirac fermion spectrum was supported by observations of the negative interlayer 
magnetoresistance. Tajima et al. reported a remarkable negative interlayer magnetoresistance under per- 
pendicular mag netic field. 20 ) The magnetic field dependence of this negative interlayer magnetoresistance 
is in good agreement with the formula given by Osada based on the zero energy Landau level. 21 ) For 
relatively large magnetic field, Osada's formula is quantitatively in good agreement with the experiments. 
In addition, in our previous publication 22 ) we reported that the interlayer magnetoresistance can be used 
to verify the tilt of the Dirac cone by extending the Osada's formula. Meanwhile, the system exhibits a 
positive magnetoresistance at low magnetic field. The origin of this positive magnetoresistance has not 
been clarified so far. 

In this paper, we argue that the positive interlayer magnetoresistance appears when the direction of the 
interlayer electron hopping is not perpendicular to the plane. Under the condition of such a non-vertical 
interlayer hopping electrons move in the plane upon interlayer hopping. Suppose electron hops from j-th 
layer to j + 1-th layer. The change of the electron coordinate in the plane implies that the n-th Landau 
level wave function in j-th layer is not orthogonal to the n'-th Landau level wave function in j + 1-th layer 
[n / n'). This inter-Landau level mixing leads to the positive magnetoresistance when \E n — E n+ \\ is less 
than the energy scale T characterizing the Landau level broadening that is satisfied at low magnetic fields. 

The rest of the paper is organized as follows. In sec. 2, we describe the model of the non-vertical 
interlayer hopping. In sec. 3, we calculate the interlayer magnetoresistance and show that there is crossover 
between a positive interlayer magnetoresistance to a negative interlayer magnetoresistance. Finally we 
summarize the result in sec.4. 

2. Model 

Following Ref.21 we calculate the interlayer conductivity by taking the interlayer transfer integral, t c , 
as a small parameter. Cartesian coordinate axes are defined in Fig. 1(a). The direction of the interlayer 
tunneling is defined in Fig. 1(b). From the crystal structure data in Ref.9 we find that the angle between 
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the interlayer tunneling direction and the z-axis is about 28 degree. This angle may depend on applied 
pressure. However, the precise value of the angle is not necessary. The value of the interlayer conductivity 
is insensitive to changes of the angle as we shall see later. The interlayer tunneling Hamiltonian is written 
as 

ifc-^/*/*£*j>,,)W*.» + *> + *- ff> 

Here the operator ipj a (x,y)^ is the creation operator for an electron with spin a at point (x, y) in the 
j-th plane. The interlayer tunneling direction, which is not parallel to the c-axis as shown in Fig. 1(b), 
is denoted by d hereafter. There is a slight translation in the x-axis in the interlayer tunneling direction 
defined in Fig. 1 (b) . We neglect that translation so that the interlayer tunneling direction is parameterized 
by c y . 




Fig. 1. (Color online) (a) The xy-coordinate system denned in the conducting planes. The a-axis is taken to be the 
x-axis. The BEDT-TTF molecules form stacks parallel to the a-axis. The z-axis is taken to be perpendicular 
to the plane, (b) The direction of the interlayer tunneling viewed from the positive a-axis. The dashed arrow 
shows the direction of the c-axis. The solid arrow shows the direction of the interlayer tunneling, a^ei, + a c e Cl 
where at and a c are the lattice constants and and e c are the unit vectors in the b and c axes, respectively. 



The interlayer current operator is found by making use of the continuity equation. We calculate the 
time derivative of the density operator under the interlayer tunneling Hamiltonian eq.(2). The result is 
equal to the divergence of the interlayer current operator in the discretize form with the minus sign. The 
effect of the magnetic field, B, is included by the Pierls substitution. We focus on the magnetic field 
perpendicular to the plane. Thus, the current operator is 

J c > = ia z t d ^2 V) ex P 

No- 
where k a = —id a {a = x,y) and a z is the interlayer distance. Hereafter we employ units where H= 1. 
Within the conducting layers the dynamics of the electrons is described by an extended Hubbard 



k y + eBx 



(3) 
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model. 6 ~ 8 ) Kobayashi el al. introduced mean fields for charge disproportionation. Since there are four 
BEDT-TTF molecules in the unit cell, four charge density mean fields were introduced. Using the self- 
consistent mean field values they found the zero gap state. In the vicinity of the zero gap point in the 
Brillouin zone, called Dirac points, the electrons are described by the effective Hamiltonian that has the 
form of a massless Dirac fermion Hamiltonian, or, more precisely, a tilted Weyl Hamiltonian. 8 ) There are 
two Dirac points in the Brillouin zone. If there are impurity potentials with extremely short interaction 
range like a lattice vacancy, there is scattering between two Dirac points. 23 ) However, it is unlikely that 
there are lattice vacancies in the BEDT-TTF molecule planes. We neglect intervalley scattering and we 
assume that two Dirac cones are degenerate. In this case we may focus on a single Dirac point. The 
Hamiltonian describing Dirac fermions in the vicinity of the Dirac cone is 



H = 

3,v 



(k x + iky^j vrjky 

In order to focus on the crossover phenomenon from the negative to positive interlayer magnetoresistance 
we have taken a simple form compared to the general form in Ref.8. In a-(BEDT-TTF)2 I3, the Dirac 
cone is tilted. 8 ) The parameter 77 is introduced in eq.(4) to describe the tilting effect. 

3. Unified Formula for Interlayer Magnetoresistance 

Now we calculate the interlayer conductivity using the Kubo formula, 

2 



i \- f (E a ) - f (Ep) 



(a\ J c 



E a — Er E a -E« + i6 



a, 1 



I E J°l dE (-JI) I («i ^ f 6 ( E - E -) 6 ( E - e p)- 



(5) 

where S is the area of the conducting layers, / = f(E) denotes the Fermi distribution function, S is a 
positive infinitesimal number, and a and /3 represent the single-body quantum states. The eigen-energy of 
the quantum state a is represented by E a . The current operator is defined by eq.(3). In order to include 
impurity scattering effects in eq.(5), we replace the delta functions, which are associated with the density 
of states, with the Lorentzian function with the half value width of F. We take T as a constant parameter 
for simplicity. 

In a magnetic field B parallel to the z-axis, the electron motion in the plane is quantized to the Landau 
levels. The derivation of the Landau level wave functions is presented in Appendix. If we take the Landau 
gauge, A = (0, Bx, 0), the Hamiltonian does not depend on y. So the wave functions take the plane 
wave form in the y-direction. Denoting the wave number by k, the quantum states are represented by 
l a ) = |i) n ) k, a) using the layer index j, the Landau level index re, and spin a = ±. We do not include the 
Dirac point index because the energy levels are independent of that index. The degeneracy with respect 
to the Dirac points is included as an overall factor in a zz . The eigen-energy is E na = E n — hbBct, where 
jj,B is the Bohr magneton. In the Zeeman energy term we have assumed that the ^-factor is equal to two 
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because spin-orbit coupling is negligibly small. 

Now we compute the matrix elements of the current operator. Using the representation of the field 
operator in terms of the Landau level wave functions <S> n (X), defined by eq.(A-13) with £ = 1/y/eB, 

4>j,o (x, y) = -7= ex P ( ik y) ® n (f + M ) c J,n,k,a, (6) 

V V n,k 

the current operator matrix elements are given by 

(j,n, k, a\ J c > | j + l,n', k', a') = ia z t d S k ^S c ^ exp (ikc y ) 

x- j dx exp (ic y eBx) + k£j $ n / Q + ki^j 

iaztc'$k,k'ficr,cr> \ r (0,0) , / A r (l,l) 

= / = j == = Y^J + S § n i nn ) 7 n> 

y (2 - <5n,o) (2 - Sn'fi) 

-r/sgn (n') J<°;V - ??sgn (n) /^;° } ] , (7) 

= j dX exp (i^x) h H _ s (x + ^2 sgn (n) V^ 3 

x/i| n ,|-«' + Y^2 s S n ( n V2A 3 Kl) . (8) 

with h n (X) being defined by eq.(A-lO) in Appendix. The function h n (X) are the eigen- functions of the 
harmonic oscillator, "H^ = p 2 /2 + A 2 X 2 /2. For the calculation of the interlayer conductivity, we need 
J(ni,n 2 ) defined through 

| (j,m,k,a\ J c > \j + l,n 2 ,k,a) \/(a z t d ) = J(n 1 ,n 2 ). (9) 

Numerical evaluations of J(ni, n 2 ) are shown in Fig. 2 for tilted (77 = 0.9) and non-tilted (rj = 0) cases with 
c y = 9.3A One crucial observation is that J(ni, n 2 ) 7^ for n\ 7^ n 2 . By contrast, for c y = 0, J(ni,n 2 ) = 
if n\ 7^ n 2 . For the tilted case, J(ni,n 2 ) gradually vanishes by increasing n 2 . For the non-tilted case, 
J(n, n + 1) is close to that for the tilted case. However, J(n, n + m) <C n + 1) for m > 2. 

For tilted cases we were unable to carry out the integral analytically. For the non-tilted case, we rewrite 
the current operator matrix elements in terms of a summation formula, 

\n'\ — \n\ 



J(n, n') = exp(- 



1 /c„\ 2 



4 w; ' \^2£ 



£1 f IfM 2 

^ . m!(|n| -m)!(|n'| - |n| +m)! \ 2\^/ 



m=max{|n| — |n'|,0} 



(10) 

Note that J(n,n) is analytically computed as J(n,n) = exp (— c 2 /(4^ 2 )) |L n (c 2 /(2^ 2 )) | with L n the n-th 
Laguerre polynomial. For n 7^ n' we need to compute J(n, ra') numerically. Figure 3 shows the magnetic 
field dependence of the current operator matrix elements. The matrix elements increase as we increase the 
magnetic field. Note that the exponential factor in eq.(10) is almost equal to one because c y /£ <C 1 for 
B < 500T. The parameter c y may depend on applied pressure. However, the range of c y is much smaller 
than the range of the magnetic length I. Since J(n, n') is expressed as a function of c y /£, the change of c y 
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Fig. 2. (Color online) The current operator matrix element J(ni,ri2) as a function of n-i for n\ = and n\ — 1 at 
B = 3T (cy/£ = 6.3 x 1CT 2 ). Both of tilted and non-tilted Dirac cone cases are shown. 



is unimportant compared to the applied magnetic field change. 

For the calculation of the interlayer conductivity, we need to include all matrix elements between 
different Landau levels: If we include n Landau levels, we need to compute n(n — l)/2 matrix elements. 
Since the expression of the matrix element for the case of the tilted Dirac cone contains the one-dimensional 
integral, the numerical calculation is heavy because it turns out that we need n ~ 100 to observe the 
positive interlayer magnetoresistance region at low magnetic fields. In order to focus on the effect of the 
non- vanishing matrix elements we use the simpler expression eq.(10) for the non-tilted case in the following. 

The interlayer magnetoresistance for C = 20KT~ 1 / 2 and T = 10K is shown in Fig. 4 . For B > 10T, 
there is a positive magnetoresistance region. This positive magnetoresistance is associated with the Zeeman 
energy effect as discussed by Osada. 21 ) Our new finding is that the appearance of the positive magnetore- 
sistance region in very low magnetic field. For B < 0.2T, there is another positive magnetoresistance 
region. From the comparison with the case of J(n,n') — > 8 nn i (dotted line in Fig. 4), it is clear that this 
positive magnetoresistance region appears because of the presence of the non-vanishing current operator 
matrix elements between different Landau levels. 

In order to make clear how the peak position is determined, we calculated the interlayer magnetore- 
sistance for different values of T (Fig. 5). The peak appears at B p ~ (T/C) 2 . The appearance of the peak 
at B = Bp is understood as follows. For B > B p , the zero-energy Landau level, n = 0, is well separated 
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Fig. 3. (Color online) The magnetic field dependence of J(ni 1 n 2 ) for different (n 1 ,n 2 ) 



from the other Landau levels because E\ — Eq > V. Therefore, the non-vanishing J(0, 1) does not play an 
important role. For B < B p , there are large overlaps between Landau levels. The Landau levels are mixed 
up upon interlayer hopping processes because of non- vanishing J(n,n'). The reason why B p is not exactly 
equal to (T/C) 2 is that the Landau levels are not equally separated in Dirac fermion systems. That is, 
E n+ i — E n ^ Ei — Eq for n^0, though the condition for the n = Landau level and the n = 1 Landau 
level plays the dominant role in the determination of the peak position. The fact that B p < (T/C) 2 is 
consistent with this picture. 

Now we analyze the experimental data in Ref.20. The key fact is that the Landau level broadening 
arises from the temperature effect and the impurity scattering effect. (Since the peak is observed for small 
magnetic fields, we may neglect the Landau level broadening effect described by the self-consistent Born 
approximation. 24 ') At high temperatures, the Landau level broadening mainly arises from the temperature 
effect. So we may neglect the impurity scattering effect. The Landau level broadening factor V is given 
by r = ksT. Therefore, if we apply the formula B p = (T/C) 2 = (ksT/C) 2 to the experimental data, 
we expect that ksE / s/Bp approaches a temperature independent constant value. The constant value is 
approximately equal to C. In order to estimate C, we evaluated the peak position, B p from the experimental 
data (not shown), and plotted fc^T '/ ' \f~B~ p in Fig. 6. We see that this quantity approaches a temperature 
independent constant value. From this plot, we found C ~ 10KT -1 / 2 . By contrast, at low temperatures the 
Landau level broadening arises from the temperature independent impurity scattering effect. Therefore, we 
expect that 1/ 'y/Bp becomes temperature independent at low temperatures. In Fig. 6, 1/ ' \f~B~ p as a function 
of temperature is shown. We see that the data is almost temperature independent at low temperatures. 
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Fig. 4. (Color online) Interlayer magnetoresistance for C — 20KT 1 / 2 and T — 10K at zero temperature (solid 
line). The dotted line shows the case of J(n,n') — > 8 n ^ n i. The inset shows the magniheation around the peak. 
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From this analysis we evaluated r as T ~ 3K. 
4. Summary 

In this paper, we have derived a formula describing the crossover from the positive to negative mag- 
netoresistance in the layered Dirac fermion system assuming a non- vertical interlayer tunneling. Although 
the tunneling direction may depend on applied pressure, the effect of the change of the tunneling direction 
is negligible compared to the change of the magnetic field because the magnetic length is much longer 
than the change of the electron coordinate in the plane. For a non-vertical interlayer tunneling, there are 
inter-Landau level mixing upon interlayer tunneling. This inter-Landau level mixing leads to the positive 
magnetoresistance at low magnetic fields. (On the other hand, the positive magnetoresistance associated 
with the Zeeman energy is observed at high magnetic fields as pointed out in Ref.21.) The inter-Landau 
level mixing effect plays an important role when the Landau level broadening factor is less than the Landau 
level energy gap. The crossover magnetic field is approximately given by B p ~ (T/C) 2 . We have applied 
this formula to the experimental data. 20 ) The Dirac fermion spectrum parameter C is evaluated to be 
C ~ lOKT -1 / 2 . It would be interesting to compare the value of C evaluated in this paper with that 
evaluated by other analyses. 
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Appendix: Landau level wave functions 

In this appendix, we derive the Landau level wave functions for tilted Dirac cone Hamiltonian (4). Under 
the magnetic field B = (0, 0, B), taking the Landau gauge, A = (0, Bx, 0), the Schrodinger equation is 

rj (jt y + eBx^j k x — i (k y + eBx^ 
k x +i(ky + eBx^j r\ (jt y + eBx^j 

Since the Hamiltonian has the translation symmetry in the y-direction, we write the wave function ip(x, y) 
as 

^ (x, y) = — \= exp (iky) </> (x) . (A-2) 
After substituting this form, we have the one-dimensional Schrodinger equation, 

4> (x) = E(f> (x) . (A-3) 



i> (x, y) = (x, y) . (A-l) 



r)(k + eBx) -i± -i(k + eBx) 



V —i-^ + i(k + eBx) rj(k + eBx) 
Introducing X = x/l + kl with t = 1/y/eB the magnetic length, (A-2) is rewritten as 

(M + 7?x) $ (X) = e$ (X) , (A-4) 
where <E> (X) = y/£<f>(x), e = El/v, and 

( -i-fir-iX \ , s 

M = . „ .„ dX „ • (A-5) 



-i&+iX 



dX 

Now we subtract the term, rfX&(X), from the both hand sides of eq.(A-4), 

M$ (X) = (e- r]X) $ (X) . (A-6) 

This is the key step to solve the Schrodinger equation. Multiplying the operator M from the left hand 
side, and using 

M(e-r}X) = (e-r}X)M+ I ° tV j , (A-7) 

\ IT) J 



we obtain 



dX 2 



TZji£ 2 - 1 irj 



*(X)=\ ^ l ^ )$(X). (A-9) 



*V T^ e2 + 1 
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The operator in the square bracket in the left hand side has the form of the harmonic oscillator with the 



(A-10) 



frequency A = yl — rj 2 . The eigen- values are (2n + 1)A and the eigen- functions are 

with X' = X + r]e/(l — rf) and H n the n-th Hermite polynomial. Making use of this result, we solve (A-9). 
The eigen-functions with the eigen-energy e are 

1 / 1 + A ' 
—irj 



V2 + 2A 



K-l (X') 



(A-ll) 



for n > 1 and 



1 



irj 



(A-12) 



V2 + 2X \ l + A / 

Note that these degenerate eigen-functions are for eq.(A-9). In order to find the eigen-states of eq.(A-4), 
we take a linear combination of these eigen-functions. After substituting the expression into eq.(A-4) and 
then after some algebra, we find the eigen-functions of eq.(A-4), 

1 



IT] 

1 + A 



h\ n \ (X n ) - isgn(n) 



1 + A 

—irj 



h\ n \-i (X n ) 



with 



Xn 



x + 



1 — rj- 



sgn (n) y / 2A 3 



n 



(A-13) 



(A-14) 



and the eigen-energy, e„ = sgn (n) ^2A 3 jn[. 
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